Centerline and tree branch skeleton determination for virtual objects

ABSTRACT

In accordance with the present invention, a method for determining a centerline through a region of interest in a 3D image dataset is provided. The method includes identifying the boundaries of the region of interest and identifying the endpoints of the region of interest. For those points within the boundaries, a penalty value which is a function of the proximity of the point to a boundary is determined. A centerline is then identified by the path connecting the endpoints which has the minimum penalized distance wherein the penalized distance reflects the actual accumulated pathlength and the penalties associated with the points along the path. From the centerline, branches of a complete skeleton can be established by determining branch endpoints and then finding the minimum penalized distance from each endpoint the centerline or another intersecting branch.

STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH

The subject matter of this application was funded in part by theNational Institute of Health, contract number CA79180 and the U.S.Office of Naval Research, grant number N000149710402. From these grants,the U.S. government may have certain rights to the invention.

CROSS REFERENCE TO RELATED APPLICATIONS

The present application is a U.S. National Stage Entry of InternationalPatent Application No. PCT/US01/30858 filed Oct. 2, 2001, the entiredisclosure of which is incorporated herein by reference. The presentapplication also claims priority from U.S. patent application Ser. No.60/237,311 filed Oct. 2, 2000, the entire disclosure of which is alsoincorporated herein by reference.

FIELD OF THE INVENTION

The present invention relates generally to computer-based threedimensional visualization and more particularly relates to improvedmethods for virtual navigation including skeleton generation.

BACKGROUND OF THE INVENTION

In many computer-based visualization systems, it is desirable to be ableto define a skeleton or centerline through the object being viewed. Sucha skeleton provides a reference path for virtual navigation, such as invirtual colonoscopy. Similarly, accurate length measurements andnavigation through other organs, such as the aorta, also require anobject skeleton. Of course, the need to define an accurate centerline ortree is not limited to medical imaging. Other fields, such as virtualengineering and architecture would also benefit from accurate methods ofdefining a skeleton in tree like objects or generating a centerline inlinear objects in various visualization applications.

Various techniques are known for generating a center line for an objectbeing visualized. For example, in the virtual examination of organs,such as the colon, a centerline can be manually marked by theradiologist while reviewing image data, such as CT, MRI or ultrasounddata. However, because of the high level of manual intervention involvedin this technique, it is expensive and generally not preferred.

An automatic method of generating a centerline is referred to astopological thinning, or “onion peeling.” While there are numerous knownvariants of this technique, in general, volume units (voxels) from aregion of interest are removed layer by layer from the boundary untilwhat remains is a connected line of predetermined thickness, such as onevoxel. While this technique has generally proven effective atdetermining a centerline in objects such as the colon, it has thedisadvantage of being computationally expensive. In addition, withcertain geometry's, the centerline that is determined by topologicalthinning is not always the centerline which would be intuitivelydesired.

Other techniques for generating a centerline through a virtual objectinclude the use algorithms which attempt to determine the shortest paththrough an object, such as the Dijkstra algorithm, which is disclosed inthe article “A Note on Two Problems in Connexion with Graphs,” Dijkstra,Numerishe Mathemetik, vol. 1, pp 269-271, 1959, which is herebyincorporated by reference in its entirety. The Dijkstra algorithm findsthe global minimal weight path in a weighted graph with nonnegativeweights. This algorithm generally includes two phases: first, a distancefrom a source field (DSF) is created by labeling all graph vertices withthe shortest distance from a single source to those vertices. Second,the shortest path is determined by tracing back from the farthest nodeto the source node. Thus, use of this algorithm requires mapping theimage data to a graph form. An example of such a mapping is illustratedin FIG. 8 a. The Dijkstra algorithm has a draw back in the case of sharpturns occurring in the object. In such as case, the algorithm in seekingthe shortest path tends to “hug the corner” and the centerline shiftsfrom the center towards the boundary.

A skeleton of a virtual object can be defined in the context of a numberof desired properties. As used herein, a skeleton should be tree-shapedand composed of single voxel paths through the branches. Thus, theskeleton represents a 1D curve in a 3D space rather than a 2D manifold.The skeleton should remain within the boundary of the region of interestand more preferably should remain in the “center” of the shape. Forobjects of arbitrary form, the concept of center may require a balancingbetween maximizing the distance from all of the boundaries of theboundaries of the region of interest and minimizing the extent that thecenterline takes on a tortuous path. This balancing implies thedesirability of applying a form of shortest path or union of shortestpaths in deriving a skeleton.

Although a number of techniques are known for generating a skeleton ofan object, each of the known techniques have inherent drawbacks. Thus,there remains a need for improved techniques for generating centerlines,skeletons and or navigation paths within virtual objects.

The centerline is often used as the primary path for navigating througha virtual object. For example, in the case of virtual colonoscopy, theuser is generally allowed to navigate along the centerline of the colonin either direction and can navigate off the centerline to exploreregions of interest. Thus, virtual colonoscopy can provide the viewerwith a more flexible view of the inner luminal surface of the colon thancan be achieved with conventional optical colonoscopy. This isaccomplished due to the fact that virtual reality and computer graphicsallow the virtual camera to view from any arbitrary viewpoint and withany arbitrary view-direction, while optical colonoscopy is limited to apath in one direction along the colon (due to the fact that the opticalcamera is placed at the end of a probe with limited mobility). However,with the flexibility that is provided by virtual examination comes a newchallenge. Specifically, the difficulty in following orientationchanges. The colon is a very tortuous organ, and, as a result, whenfollowing along the colon in 3D users can become disoriented. This riskis especially present when the user veers off the centerline andexamines a structure close to the colon wall. Often times, whenreturning to the centerline, the user may head back in the directionthey previously came from. Since the view is from the oppositedirection, it is usually not readily apparent that the flight is headedin the reverse direction. Thus it would be desirable to provide a methodof reducing the opportunity for a viewer to become disoriented whilenavigating through the virtual object without reducing the desirableflexibility which can be achieved in the virtual environment.

SUMMARY OF THE INVENTION

It is an object of the present invention to provide an improved methodfor generating a centerline, and where applicable a skeleton, of avirtual object.

In accordance with the present invention, a method for determining acenterline through a region of interest in a 3D image dataset isprovided. The method includes identifying the boundaries of the regionof interest and identifying the endpoints of the region of interest. Forthose points within the boundaries, a penalty value which isproportional to proximity of the point to a boundary is determined. Acenterline is then identified by the longest Euclidean path connectingthe endpoints which has the minimum penalized distance wherein thepenalized distance reflects the actual accumulated pathlength and thepenalties associated with the points along the path.

The method preferably includes steps to initially identify points whichare near the centerline of the object. In this case, the steps ofdetermining a penalty value and determining a path described above canbe performed only on those points which are near the centerline. Thisreduces the dataset which is to be processed thereby improvingprocessing cost and speed.

The operation of identifying points near the centerline can beaccomplished by: computing a distance from boundary field (DBF) for eachof the points in the region of interest; determining a gradient fieldfor each point in the region of interest; identifying regions of pointshaving a non-uniform gradient; and connecting the regions having anon-uniform gradient.

In addition to a single centerline, it is also desirable to extend themethod to determine a complete skeleton of a branched object. This canbe performed by identifying the endpoints of the branches extending fromthe centerline, and for each identified endpoint, determining a branchpath connecting the endpoint to another branch which has the minimumpenalized distance. Again, the penalized distance for each branchreflects the actual accumulated pathlength of the branch and thepenalties associated with the points along the branch path.

BRIEF DESCRIPTION OF THE DRAWING

Further objects, features and advantages of the invention will becomeapparent from the following detailed description taken in conjunctionwith the accompanying figures showing illustrative embodiments of theinvention, in which

FIG. 1 is a flow chart illustrating a preferred method for generating askeleton of a virtual object in accordance with the present invention;

FIG. 2 is a graphical representation of an outline of a region ofinterest along with a numerical representation of the Euclidean distancebetween each inside voxel and the object boundary or distance fromboundary field (DBF);

FIG. 3 is a graphical representation of gradient field vectors forvoxels in a region of interest;

FIG. 4 is a graphical representation illustrating the flagging ofnonuniform gradient neighborhoods from the diagram of FIG. 3.

FIG. 5 is a flow chart illustrating a method for labeling nonuniformgradient vector field vectors;

FIG. 6 is a flow chart for connecting flagged voxels of a centerline;

FIG. 7 is a graphical representation illustrating the connection offlagged voxels along the gradient vector field vectors;

FIGS. 8 a and 8 b are a pictorial diagram illustrating a mapping of avoxel grid and neighbor relations to an undirected graph on which theDijkstra algorithm can be applied, wherein FIG. 8 a represents aconventional shortest distance source field and FIG. 8 b represents apenalized distance from the source field, which is used in the presentinvention;

FIG. 9 a illustrates the generation of a centerline through an aorta;

FIG. 9 b illustrates the generation of a complete tree skeleton throughan aorta, in accordance with the method of FIG. 1.

FIG. 10 is a cross section of a region of interest with a centerline andparameters for establishing the direction of travel along thecenterline.

Throughout the figures, the same reference numerals and characters,unless otherwise stated, are used to denote like features, elements,components or portions of the illustrated embodiments. Moreover, whilethe subject invention will now be described in detail with reference tothe figures, it is done so in connection with the illustrativeembodiments. It is intended that changes and modifications can be madeto the described embodiments without departing from the true scope andspirit of the subject invention as defined by the appended claims.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

The present invention provides an improved and computationally efficientmethod of automatically generating a centerline in a lumen or a skeletonin a more complex object. The method uses the concept of penalizeddistance fields to augment known shortest path algorithms, such as theDijkstra shortest path. The examples referred to are generally in thecontext of medical examination, such as virtual colonoscopy, but areapplicable to any application of computer based virtual navigation.

As noted above, a problem with the known Dijkstra shortest pathalgorithm is evident when the path through the object turns sharply. Insuch a case, without some form of embellishment or modification, theDijkstra shortest path will tend to hug the corner and deviate from thecenter of the object. In the present methods, this is corrected by useof the notion of a penalty field. FIG. 8A illustrates a distance fieldconnecting voxels in accordance with the conventional Dijkstraalgorithm. Referring to FIG. 8B, the present method supplements theneighbor edges of the Dijkstra algorithm by adding “penalty edges” whichare weighted to provide an increased weight or penalty for those voxelsnear the object boundary. In the representation of FIG. 8B, there are 27vertices per voxel: one central vertex and 26 penalty vertices, thateach share a penalty edge with the central vertex. The penalty edgeshave a weight equal to one half the penalty associated with includingthat voxel in the path. In this case, neighbor edges now connect topenalty vertices. This results in a graph that is a singly connectedcomponent where edges have positive weights to enable the Dijkstraalgorithm to find the globally minimal shortest path. The costassociated with a path is the accumulated piecewise Euclidean distanceof the path plus the sum of all of the penalties associated with thepath. The centerline is defined as the minimum cost path through thepenalized distance field.

FIG. 1 is a flow chart depicting an method for determining a skeleton ofa virtual object. The method assumes that image data of an object hasbeen acquired, such as CT or MRI data, and that this data hastransformed into the 3D volume domain and has been segmented. Acquiring2D image data, such as CT, MRI, ultrasound, and the like andtransforming the 2D data into a 3D representation having a plurality ofvoxels is well known in the art, as are various techniques forsegmenting such 3D data. The data is thus presented as a 3D rectilineargrid which can be referred to as the volume wherein each volumetricsample point is called a voxel.

For the purposes of determining the centerline, the segmented image datacan be subjected to a binary mask which delineates the object boundaryand labels those voxels within the volume of the boundary as “inside”voxels. (step 102). The segmentation techniques disclosed in connectionwith electronic colon cleansing in U.S. patent application Ser. No.09/777,120, filed on Feb. 2, 2001 and entitled “System and Method forPerforming a Three-Dimensional Virtual Examination of Objects, Such asInternal Organs,” which is hereby incorporated by reference, can be usedin this regard. With the boundary defined by the binary mask, the regionof interest can be cropped using a bounding box which encloses thosevoxels tagged as “inside.” (Step 104). While not essential to thepractice of the invention, by cropping the data in step 104, asignificant reduction in the data set is achieved and the subsequentsteps are rendered more efficient. For example, in the case of medicalimaging scans, the volume size is reduced on the order of 30-50%.

Following cropping, the Euclidean distance between each voxel flagged as“inside” the volume and the boundary is determined and recorded. (Step106). As illustrated in FIG. 2, this results in a mapping of thedistance from boundary field (DBF). One method of determining theEuclidean distance is a four pass algorithm disclosed in the article“New Algorithms for Euclidean Distance Transformation of ann-Dimensional Digitized Picture with Applications,” T. Saito et al.,Pattern Recognition, vol. 27, no. 11, pp. 1551-1565. However, otheralgorithms for calculating the distance from the respective voxels tothe boundaries can also be employed.

Preferably, for each voxel within the boundary of the region ofinterest, a central difference gradient is calculated to define agradient vector field (GVF) of the distance from boundary fields (DBF).(Step 108). This can be achieved by comparing each voxel with only sixneighboring voxels. The GVF at point (x,y,z) is found byG_(xyz)=(G_(x),G_(y),G_(z)) where G_(x)=V_(x+1,y,z)−V_(x−1,y,z) andG_(y)=V_(x,y+1,z)−V_(x,y−1,z) and G_(z)=V_(x,y,z+1)−V_(x,y,z−1). Highernumbers of neighboring voxels can be used in this calculation, such asby using a 26-voxel neighborhood Sobel filter, however increasing thenumber of voxels tends to slow the calculation without improving theresulting centerline. The graphical representation of FIG. 3 illustratesa GVF for the representation of the DBF of FIG. 2. In FIG. 3, eachgradient vector field vector associated with a respective voxel isillustrated as an arrow 300 with its origin at the voxel position.

After the GVF is established, the gradient vector field is analyzed todetermine the characteristics of local neighborhoods within the gradientfield vectors. The neighborhoods can be defined as 2×2×2 voxel cells.Referring to FIG. 4, the neighborhoods can be characterized in sixclasses: local maxima near the centerline 402, local minima near thecenterline 404, uniform areas near the centerline 410, uniform areas offthe centerline 412, local maxima off the centerline 406 and local minimaoff the centerline 408.

Local maxima on or near the centerline 402 are defined as those regionsnear the centerline of the object in which all of the gradient vectorfield vectors of the 8 voxel neighborhood point in a direction towardsthe centerline. Thus, as illustrated by local maxima 414 and 416, all ofthe gradient vector field vectors are pointing in a different, butconverging direction.

A local minima on the centerline 404 is defined as a neighborhood whereall of the gradient vector field vectors surrounding such a minima aredirected away from this point and towards a neighboring maxima. Thevectors will then be in at least two groups of vectors pointing inopposite directions. In the case of a centerline through a lumen, eachminima is bounded by a maxima on either side.

A local maxima located off of the centerline 406 generally has the samecharacteristics as a local maxima which is located on or near thecenterline 402. The only difference being that such a local maxima ispositioned remote from the centerline. This can occur in branches ordistortions in the object off of the main centerline, such asillustrated by local maxima 418. Similarly, local minima off thecenterline, such as local minima 420 are indistinguishable from localminima on the centerline, such as local minima 422, except for theirposition with respect to the centerline.

In addition to local maximas and local minimas there are also regions ofuniform gradient vector field vectors, both on the centerline 410 andoff the centerline 424. A uniform GVF is defined as one where the GVFvectors of the neighborhood, such as an eight voxel cell, have adirection which varies less than ninety degrees from the averagedirection. In those regions which are near the object boundary, thevectors in a local neighborhood generally point in the same directionwhich is perpendicular to the boundary. There can also be regions alongthe centerline where the object widens quickly. As a result, the GVFvectors neighboring the centerline in such regions will tend to pointuniformly towards the center of this region. Thus, the method oflabeling all nonuniform gradient vector field regions does not capturethis portion of the centerline.

Given that a uniform GVF is defined as one where the GVF vectors of theneighborhood, such as an eight voxel cube or cell, have a directionwhich varies less than ninety degrees from the average direction, theprocess set forth in FIG. 5 can be used to label local uniform andnonuniform GVF vectors. In step 505 the normalized average gradientfield vector is computed for each defined neighborhood, such as 2×2×2cells of voxels, for each voxel position. Next, the dot-product betweeneach of the 2×2×2 GVF vectors and the associated average vector isdetermined and tested. (step 510). A positive dot product indicates auniform GVF, whereas a zero or negative dot-product is indicative of anonuniform GVF. Thus, to identify the nonuniform neighborhoods, thosecells having a nonpositive dot product are labeled. (Step 515).

The use of a small number of voxels, i.e., 8, in each voxel positionneighborhood allows for rapid execution time and results in the flaggingof only a small percentage of voxels inside the boundary as beingpotentially near the centerline.

The result of the flagging process is a number of disconnected groups offlagged voxels. The next step is to connect the flagged voxels along thecenterline (FIG. 1, Step 110). An overview of a procedure for connectingthe flagged voxels of the nonuniform GVF's is illustrated in FIGS. 6 and7. First, a flagged voxel position, such as local minima 705, isselected as a starting point. (Step 605). Each vector of the flaggedvoxels has a component which points along the centerline and a componentwhich points towards the centerline. Thus, a voxel by voxel path istraversed from the first flagged voxel, along the direction of the GVFvector (step 610), until a next flagged voxel is reached (step 615).Referring to FIG. 7, beginning with flagged voxel 710, the GVF vector isdirected to voxel 712. This process is repeated until flagged voxel 714of local maxima 716 is reached. Because the GVF vectors tend to pointtowards the centerline, the process of traversing along the voxels is aself correcting process. For example, voxel 710 is to the right of thecenterline and the next voxel selected 712 is to the left of thecenterline. Because the GVF vector of 712 is directed towards thecenterline, the next voxel selected will at least be closer to thecenterline or, perhaps, to the right of the centerline. The process ofsteps 605 through 615 is repeated for each group of flagged nonuniformneighborhoods.

Returning to FIG. 1, a starting point voxel within the object isselected and the farthest voxel(s) from that voxel are determined bycalculating the distance from the source field (DSF) (step 114). Suchfarthest voxels from any arbitrary voxel will represent the endpoints ofthe centerline and can be selected as a root voxel for the centerline orskeleton. The choice of starting point voxel can be improved by domainknowledge. In the case of colonoscopy, the rectum is the optimum pointto start the examination, so we wish to select a point near the rectum.An estimation of the location of this point can be determined by findingthe nearest “inside” point closest to the center-bottom of the dataset.This can be found by executing a search through the voxels near thecenter-bottom dataset voxel in shells of increasing radii until a voxeldetermined to be “inside” the colon is found.

In step 116, the penalized distance from the root voxel field (PDRF) isdetermined. Calculating the PDRF involves assigning a penalty to eachvoxel in the proposed path. The penalty, p at a voxel v is assignedbased on the DBF value at that voxel calculated in step 106 as well as aglobal upper bound of all DBF values, M (M>max(DBF)). The penalty, p canbe calculated as:p(v)=5000·[1−DBF(v)/M] ¹⁶

The term [1−DBF(v)/M] is always in the range of [0, 1] with the maximalvalues occurring for voxels near the boundary of the object. The term5000, which is chosen heuristically, is selected to be sufficientlylarge to ensure that the resulting penalty is sufficiently large toovercome any distance advantages of a straight path. The value of 5000allows skeleton segments to be up to 3000 voxels long without exceedingfloating point precision. Other functions with similar characteristicscan be chosen (e.g., small penalty towards the center and large at theboundary).

After assigning the penalty to each of the flagged voxels, the minimumcost path can be determined. (Step 118). The accumulated cost at eachvoxel along a path can be expressed as:cost(v _(k))=cost(v _(k)−1)+distance(v _(k) , v _(k−1))+penalty (v _(k))In addition the accumulated cost, the accumulated distance at a givenvoxel, D_(a)(v), can be expressed as:D _(a)(v _(k))=D _(a)(v _(k−1))+distance(v _(k) , v _(k−1))This accumulated distance can be used to identify the farthest voxel inthe PDRF, which represents the opposite end of the centerline.

From the farthest voxel, the minimum cost path back to the source voxelis determined. This is performed, in a similar manner to the Dijkstraalgorithm, by tracing the discrete centerline from the end voxel back tothe root voxel. Because of the inclusion of higher penalties for voxelsnear the boundary, the minimum cost path tends to be pushed into thecenter of the object. This minimum cost path is the discrete centerlineof the object. In the case of a tree like structure, this discretecenterline is the main branch of the discrete skeleton tree.

The discrete centerline is made up of numerous connected short straightline segments from voxel to voxel. Such a discrete “stair step”representation is not the most desirable centerline in manyapplications, such as virtual navigation. In such a case, the stair stepcan manifest itself as jitter in the virtual camera. Therefore, in manycases it is desirable to smooth the discrete centerline. (Step 120). Onepreferred method of smoothing the curve is to use an approximatingspline algorithm with adaptive error tolerance. In narrow regions of theobject, the error value should be small whereas in wider regions, moredivergence from the discrete centerline can be tolerated. This errortolerance can be expressed as a percentage of the distance from theboundary. A percentage of less than 50% insures that the centerlineremains closer to the centerline than to the boundary, thus it isguaranteed to be within the boundary. A well known B-spline curve whichinterpolates the first and last voxel and approximates the ones inbetween can be used with control points placed close to centerlinevoxels at nonuniform intervals. A number of heuristics can be used tominimize the number of control points required to acheive a particularaccuracy. For example, an error tolerance of 35% may require 17% of thediscrete centerline voxels as control points whereas an error toleranceof 50% would only require 13% of the discrete centerline voxels ascontrol points.

In the case of a lumen-shaped object, such as a colon, where only asingle centerline is desired, the smoothing operation would complete thecenterline generation process. In the case of more complex, branchedstructures such as the lungs and aorta this initial centerline willserve as the main branch in the complete tree structure skeleton of theobject.

The determination of the rest of the branches of the skeleton for acomplex object are determined in steps 122 through 130 which arerepeated for each branch of the skeleton.

The process of extending the centerline to a complete skeletondefinition starts with the concept of “rolling an adaptive sphere” downthe skeleton voxels. This, as explained below, will result in thelabeling of those voxels which are near the centerline of the branchesin the skeleton. (Step 130). The adaptive sphere is defined at eachvoxel along the skeleton in terms of its proximity to a boundary, a userapplied scaling factor and a user defined constant. This is expressedmathematically as:radius (v)=DistanceFromBoundary (v)·scale+constantThose voxels within the radius of the adaptive sphere are labeled asprocessed voxels. The combination of the terms scale and const determinethe minimum feature size that will be used to determine a skeletonbranch. For example, in a lumenal structure, such as the colon, a scalevalue of scale=3 and const=50 results in the initial centerline labelingall colon voxels and thus remains the only centerline through theobject. For branched objects, such as the aorta, the values scale=1.1and const=10 have been found to result in the identification of allsignificant vessels as branches of the skeleton.

To identify a branch along the centerline, a new penalized distance fromroot field (PDRF) is calculated for each of the voxels labeled asprocessed voxels in step 122. (Step 124). The PDRF is performed in asimilar fashion as described above in connection with step 116. Thevoxel identified as the farthest voxel from the PDRF calculation is usedas the tip of the next major branch.

From the tip voxel of the current branch, a path is determined along thePDRF for the current branch until an existing skeleton voxel is reached,i.e. the branches merge. (Step 126). As set forth above with respect tostep 118, this involves identifying the minimum cost path through thePDRF of the branch voxels.

After the discrete minimum cost path is determined for the currentbranch, the branch centerline can be smoothed. (Step 128). The smoothingoperation can be performed using the process described in connectionwith step 120 or any number of other curve smoothing techniques whichare known. The one additional constraint that should be applied is thatthe smoothing operation should not result in a discontinuity at thepoint where two branches of the skeleton merge.

After the centerline for the branch is smoothed, those voxels of thebranch centerline are added to the dataset of skeleton branches. (Step130). Steps 122-130 are repeated until the entire segmented volume ofthe object of interest has been analyzed.

FIG. 9A illustrates the results performing steps 102 through 120 on aset of image data of an aorta. The aorta is rendered substantiallytransparent and the centerline 902 of the major branch of the aorta ishighlighted. In FIG. 9B, steps 120 through 130 have been iterativelyrepeated to extend centerline 902 into a full tree skeleton of the aortain which all branch vessels have a centerline therethrough.

A further, optional method to improve the determination of bothendpoints for any segment can be employed. The method is to compute thePDF from the approximate start point and find the maximum distance pointfrom the approximate start point. The maximum distance point is thenconsidered the true end point. Another PDF is computed from the true endpoint and the maximum distance point is then considered the true startpoint.

There may be more than one connected component (or segment)automatically identified as part of the colon. Some or all of thesesegments may actually be part of the colon. If more than two segmentsexists, the second segment may be selected by finding the nearest newsegment voxel to the end of the first segment centerline. This point canbe used as the starting point for the new PDRF in the new segment.

When two or more segments have been automatically identified by thesystem and the viewer disagrees with the assigned order, the user can begiven the opportunity to select some or all of the segments (becausesome may be pockets of air in other organs not of interest). Theselection can be made by clicking on each segment in the desired orderof traversal. For example, in colonoscopy, the doctor usually clicks onthe segment with the rectum first and the segment with the cecum last.Furthermore the user is able to reverse the preferred direction oftravel within each segment so that the connection of all centerlinesappears to be one continuous path. Furthermore, the selection ofsegments in order and the selection of directions within each segmentcan be accomplished by a single selection at one segment end or theother. The system can automatically determine which end the userselected and reverse the centerline direction if necessary.

In order to facillitate guided navigation throughout a skeleton, a usercan define a start point and an endpoint for the navigation by selectingthese points on the skeleton using a graphical user interface. If thestart or end points are near the endpoints of the centerline or branch,selection need not be precise as the selected point can be snapped tothe endpoint of the branch tips. If desired, new navigation fields canbe calculated for this newly defined path and a new centerline can begenerated between the selected points to reduce any transistion atbranch junctions.

After a centerline of skeleton is determined through a virtual object,such a centerline is often used as a primary path for navigating throughthe object. As noted above, the user is generally provided with theopportunity to navigate in both directions along the centerline and canalso navigate off of the centerline to explore regions of interest.However, with this flexibility in navigation comes the opportunity forthe user to become disoriented.

To minimize the opportunity for disorientation while navigating in aninteractive mode, the present method can employ a technique to color thecenterline to indicate the direction of travel and to mark areas whichhave already been explored. The centerline is colored according towhether the user is traveling along the same direction as the currentexamination direction (i.e., rectum→cecum or cecum→rectum). Asillustrated in FIG. 10, the normal color (such as green signifies thatthe viewpoint is within 135 degrees of the closest centerline direction.The opposite color (such as red) signifies that the viewpoint has comewithin 45 degrees of the opposite-direction of the closest centerlinedirection. The determination of this coloring is preferably filtered bya predetermined threshold value to add a level of hysterisis so that thecenterline does not flash between the two colors (FIG. 1). The anglebetween the current view direction and the closest centerline directionis simply computed as the arccos of the dot product of the twodirections. To be able to easily determine the closest centerlineposition during navigation, all valid positions for the virtualviewpoint store a pointer to the closest centerline point which ispreferably defined by discrete points rather than as a continuousrepresentation for quick lookup.

The PDF has another important use other than the creation of acenterline. The PDF can be used to translate the viewpoint of the userthrough each segment toward each segment centerline endpoint. The PDFfield has the property that at every point the gradient simultaneouslypoints away from the boundary and towards the segment centerlineendpoint. By simply translating the viewpoint along the gradientdirection until the viewpoint reaches the global minimum value, theviewpoint is smoothly drawn towards and along the centerline. This makesit easy for the user to efficiently traverse the colon without gettinglost. Furthermore, the view direction can be rotated toward thedirection of travel so that the user is looking where he is going. Tosmooth the transition from the current view direction to the new travelview direction, the view direction is combined with the new travel viewdirection with a weighted vector addition at each time step. Besidesautomatic navigation, the user can be permitted to freely navigatewithin the boundary of the object away from the centerline. The value ofthe DFB provides a measure of the virtual camera's proximity to theobject boundary. If the viewpoint is within some threshold distance fromthe boundary, then a force can be activated which smoothly pushes theviewpoint “away from” the boundary. The “away from” direction can beestimated by using the gradient of the DFB. This gradient yields themost likely direction which points away from the boundary. The magnitudeof the force can be determined by a function which is greatest closestto the boundary and decreases away from the boundary.

The cost field computation which employs penalties which areproportional to the proximity to the boundary of an object discussedabove with respect to centerline generation can also be used in guidednavigation through an object. Thus, the cost of a chosen navigation pathcan be the accumulated Euclidean distance of the chosen path plus theaccumulated penalty cost of the path. Regardless of the starting pointfor a guided navigation, the penalty fields will tend to push thenavigation towards the centerline and then direct the navigationsmoothly down the centerline towards the selected seed point, which istypically chosen as an end point of the object.

The present methods can be performed on any number of computerprocessing platforms. It is expected that the increased efficiency thatis provided by these techniques make them well suited for use onpersonal computer systems. For example, an IBM compatible Wintel PCoperating under the Windows2000 operating system and having a single 1GHZ CPU resulted in very reasonable processing times, such as generatinga colon centerline less than two minutes, and the aorta skeleton of FIG.9B in 36 seconds.

Although the present invention has been described in connection withspecific exemplary embodiments, it should be understood that variouschanges, substitutions and alterations can be made to the disclosedembodiments without departing from the spirit and scope of the inventionas set forth in the appended claims.

1. A computer-based method for determining a centerline through a regionof interest in a 3D image dataset comprising: identifying the boundariesof the region of interest in the 3D image dataset; identifying theendpoints of the region of interest; determining a gradient field foreach point in the region of interest; identifying regions of points inthe region of interest where the gradient field exhibits a non-uniformgradient; and connecting the regions having a non-uniform gradient toidentify a plurality of points proximate to the centerline of the regionof interest; for those plurality of points proximate to the centerline,determining a penalty value which is a function of proximity of thepoint to a boundary; and determining a path connecting the endpointswhich has the minimum penalized distance wherein the penalized distancereflects the actual accumulated pathlength and the penalty valuesassociated with the points along the path.
 2. The method of claim 1,further comprising smoothing the path connecting the endpoints.
 3. Themethod of claim 1, wherein at least one of the endpoints is selectedbased on prior knowledge of the object.
 4. The method of claim 1,further comprising: identifying the endpoints of branches from thecenterline; for each identified endpoint, determining a branch pathconnecting the endpoint to another branch which has the minimumpenalized distance wherein the penalized distance reflects the actualaccumulated pathlength and the penalties associated with the pointsalong the branch path.
 5. The method of claim 4, further comprisingcomputing a distance from boundary field for each of the points in thebranch paths.
 6. The method of claim 4, further comprising, for eachbranch, identifying points which are near the centerline of the branchand performing the steps of determining a penalty value and determininga path only on the points near the centerline.
 7. The method of claim 6,wherein the operation of identifying point near the centerline of thebranch further comprises: determining a gradient field for each point inthe region of interest; identifying regions of points having anon-uniform gradient; and connecting the regions having a non-uniformgradient.
 8. The method of claim 7, further comprising smoothing thecenterline of each branch.
 9. A method of generating a centerlinethrough a virtual object comprising: acquiring a 3D image dataset of thevirtual object, the image dataset comprising a plurality of voxels;identifying the boundaries of a region of interest in the 3D imagedataset; computing a distance from boundary field for each of the voxelsin the region of interest; determining a gradient voxel field for eachvoxel in the region of interest; identifying regions of voxels having anon-uniform gradient; connecting the regions of non-uniform gradient;from the connected regions, identify a root voxel; for the connectedvoxels, calculate a penalized distance from the root voxel; anddetermine the minimum penalized distance from the root voxel to thefarthest voxel from the root voxel.
 10. The method of claim 9, furthercomprising smoothing the path defined by the minimum penalized distance.11. The method of claim 9, further comprising: determining the endpointsof branches; for each branch endpoint, determining a branch pathconnecting the endpoint to another branch which has the minimumpenalized distance wherein the penalized distance reflects the actualaccumulated pathlength and the penalties associated with the pointsalong the branch path.